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Abstract 

It is shown that the algorithm introduced in [T| and conceived to 
deal with continuous degrees of freedom models is well suited to com- 
pute the density of states in models with a discrete energy spectrum 
too. The q = 10 D = 2 Potts model is considered as a test case, and it 
is shown that using the Maxwell construction the interface free energy 
can be obtained, in the thermodynamic limit, with a good degree of 
accuracy. 

1 Introduction 

In recent years Wang-Landau sampling (WLS) [21 [3] has became a standard 
tool in the area of Monte Carlo investigation of statistical systems and has 
been successfully applied to an ever increasing number of models and situa- 
tions. Since its introduction, a lot of studies have been carried over in order 
to investigate its convergence properties and to refine it [U E]. 

The WLS finds its natural application to discrete energy models and its 
use for the study of continuous energy models deserves some care, since it has 
been found, for example, that a naive "bin-discretization" of a continuous 
energy leads to difficulties that need to be handled with care [6j [TJ E] • 

In a recent publication [1J a variation on the WLS theme, conceived for 
models with a continuous energy spectrum and based on a local linearization 
of the logarithm of the density of states (DOS), has been proposed. The 
purpose of this paper is to demonstrate, at numerical level, that the algorithm 
proposed in pQ, which will be referred to as the LLR sampling, or LLRS, is 
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effective for models with a discrete energy spectrum too. There are in fact 
no a priori reasons the algorithm proposed in [TJ should not work for discrete 
energy models, at least in the infinite volume (or thermodynamic) limit, in 
which the energy density, and accordingly the DOS, becomes in any case 
continuous. Clearly, on general ground, at finite volume discrepancies with 
respect to a standard WLS, at the level of 0(1/ L a ), where a depends on the 
specific observable and L is the linear size of the lattice, are expected, but it 
will be shown that for L large enough these finite size effects are small and 
do not affect a safe extrapolation of physical quantities to the infinite volume 
limit. 

The motivation for applying the LLRS to discrete energy models stems 
from the fact that in some cases the derivative of the logarithm of the DOS 
respect to the energy (i.e. the inverse of the microcanonical equilibrium 
temperature at fixed energy) is all is needed: in such cases the LLRS gives the 
answer in a direct way and the overall computational cost can be dramatically 
decreased respect to WLS. 

In the following the algorithm described in pQ will be introduced. Then 
as a first test the Ising model in D = 2 will be addressed and the algorithm 
will be used to numerically extrapolate the microcanonical equilibrium tem- 
perature for a given energy to the thermodynamic limit; the comparison with 
the analytically known result will show that the LLRS does not produce sys- 
tematic errors in the infinite volume limit. The third step will be the study 
of the interface free energy in the D = 2, q = 10 Potts model by means of 
the so called Maxwell construction, with the help of the WLS. Afterwards 
the same computation will be repeated with the LLRS, using larger lattices. 
At the end some conclusions will be drawn. 

2 The algorithm 

Generally speaking, the canonical partition function of a statistical system 
can be expressed as an integral (or a sum, for a discrete energy model) over 
all allowed energies: 



where (3 = 1/T is the inverse temperature and g(E) is the DOS (the Boltz- 
mann normalisation factor is set to 1 through all this paper). It is observed 
that in general, maybe neglecting the boundaries of the domain of defini- 
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tion, the logarithm of the DOS (i.e. the microcanonical entropy) is a smooth 
function of the energy (equivalently of the energy density); so smooth, in 
fact, that becomes meaningful to locally apply a linear approximation and 
to write 

log g(E) = a(E )[E - E ] + c{E ). (2) 

It is understood that this approximation works only if E is in the neighbour- 
hood of E ; to introduce the notation it can be written that eq. works 
well for E - 5E/2 < E < E + 5E/2, with a suitable choose of 8E. 

The general strategy outlined in [I] will be now briefly recalled. Assuming 
the energy is restricted to the interval Eq — 5E/2 < E < E + SE/2, the 
truncated expectation value of a general observable f(E) can be defined: 

1 rE +SE/2 

((f(E Q )))(a) = - dEf(E)g(E)e- aE , (3) 

Al Je -5e/2 

where the normalisation M is given by 

M = / dEg(E)e~ aE . (4) 

JE -SE/2 

If eq. 02]) was a perfect approximation then a in eq. ([3]) could be chosen in 
order to cancel exactly g{E); in this flat distribution would remain, 

so that 

((E))(a) = E for a = a(E ). (5) 

Assuming now that a good guess for a(E ) (call it a n ) is known, ignoring 
higher order corrections and defining E* = E — E , the truncated expectation 
value of E* is readily obtained: 

((E*))(a n ) = ^f[a(E ) - a n ] + 0(x 3 5E), (6) 

where x = [a(E ) — a n ]5E. Eq fl6]) can now be solved for a(E ) in order to 
obtain a better approximation a n+ i\ 

12 

a n+1 = a n + ——^((E*))(a n ); (7) 

the game can in principle be iterated up to convergence. 

Obtaining the truncated expectation value ((E*))(a n ) is quite easy: all is 
needed is an entropic sampling [10], restricted to the relevant energy range, 
with g(E) given by eq. (T2J) and a(E ) replaced by a n . 
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As a matter of fact, as it is noted in PQ, it is quite unuseful to obtain a 
really precise estimate of {(E*) )(a n ) in order to increase the precision on a n+ i. 
In any case the value of a n , after a suitable thermalisation period, will start 
to fluctuate around its mean value, but the amplitude of the fluctuations (i.e. 
the variance of a) is of course a function of the number of entropic sampling 
sweeps used to compute ((E*))(a n ). A good compromise consists in using a 
reasonable number of entropic sampling sweeps N e , and to average a n over 
1 < n < N s after N t thermalisation steps. 

The above line of reasoning is enough to define the following algorithm 
to compute 9 log g(E)/dE for a given value of Eq and a given choice of 5E: 

1. guess a reasonable starting value for a n (n = 0); 

2. compute {(E*))(a n ) performing N e sweeps of entropic sampling with 
log g(E) = a n E (the sampling is restricted to the range [E —SE/2, E + 
5E/2]; configurations whose energy is out of this range are simply not 
taken into account in the averaging process); 

3. compute a n+1 by using eq. (J7J); 

4. n i— n + 1 and iterate; 

5. discard the first N t values of a n and average the remaining N s values 
to obtain an estimate of a(Eo), which can be directly interpreted as 



The Ising model in 2D is the first-choice testbed because of the existence of 
an analytic solution against which numerical results can be checked. In this 
case, however, the analytic result at finite volume pJJ] are of little help, since 
the LLR sampling as defined above will show finite size effects respect to the 
"true" finite size results for g(E). This is by no means a problem, since the 
check will be done against results in the thermodynamic limit. 
The hamiltonian of the Ising model is taken to be 



(d\ogg(E)/dE) E=Eo . 
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Ising model 



ft = o El 1 -^^i] 



(8) 



where the sum is over nearest-neighbours sites and cr^ = ±1. 
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The LLR sampling denned in the previous section was used in order 
to compute a(E ), with E = 0.5 V, at several lattice size L, in order to 
proceed to the infinite volume limit. In this limit a(Eo) should become 
(dS(E) I dE) e=e , a quantity for which an analytic answer is of course known. 
For each value of L the following procedure has been followed: 



L 


■^run 


a(E/V = 0.5) 


64 


50 


0.754679(7)* 


72 


50 


0.755187(6)* 


80 


37 


0.755549(7)* 


96 


50 


0.756020(5) 


112 


40 


0.756299(4) 


128 


40 


0.756494(3) 


160 


20 


0.756716(4) 


192 


20 


0.756837(5) 


240 


30 


0.756933(2) 


320 


18 


0.757005(2) 


extrap. 




0.757105(2) 



Table 1: Simulations parameters and results for Ising model. Values marked 
with an * have not been used in the extrapolation to the infinite volume limit. 

• after N t = 100 thermalisation steps, a has been averaged over N s = 
1000 measures, using 5E = L; 

• each measure, see eq. ([7]), has been obtained with N e = 1000 entropic 
sampling sweeps, where one entropic sampling sweep consists in trying 
for V = L 2 times to flip a randomly chosen spin; 

• The procedure has been repeated N run times in order to compute errors. 
Lattice sizes, N run and results are shown in table [TJ 

A few comments about the parameters used in the procedure delineated 
above are in order: 

• Starting with ao = 0.8, the convergence to the equilibrium value of a 
is really fast and N t = 100 is more than enough for all lattice sizes 
considered; 
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• N e is in any case taken proportional to the volume of the system in 
order to ensure that in the entropic sampling process each spin of the 
system has a chance to be flipped many times; it has to be taken into 
account that the "flipping acceptance" (in the present case in which 
E/V = 0.5), i.e. the probability that the flip of a given spin during 
the entropic sampling process is accepted, is about 0.37; 

• SE is taken proportional to L in order to ensure that in the infinite 
volume limit this range do not shrink to zero when measured in "phys- 
ical" units; the relevant point is that 5E/V will go to as 1/L. Taking 
SE = L seems to be the most simple choice, and with this choice the 
linearization (J2J) is satisfied to an high degree of accuracy. 

Moreover it has to be noted that a discrete energy model requires sums to 
be written instead of integrals in eqs. (JSJ) and (j3J); accordingly, for the Ising 
model eq. is modified as follows: 

an+1 = an+ A(SEH(SE/ {E * ))M (9) 

Several lattices, ranging in size from L = 64 up to L = 320 have been 
considered, and in figure [T] the thermodynamic extrapolation in 1/V (not 
considering the three smallest lattices, but their inclusion hardly changes the 
result) is shown. The final result, 

a(E ) = 0.757105(2), (10) 

has to be compared with the analytic answer, for which the convention de- 
fined in (JSJ) gives, up to six digits, the value 0.757106. The almost perfect 
agreement shows that in the thermodynamic limit there are no systematic 
effects, at least in the high statistical precision reached. The extrapolation 
shows that the finite size effect on (dS(E)/dE) for a given value of E is 
proportional to 1/L 2 . The whole computation took about three days on a 
desktop computer (quadcore). 



4 The Maxwell construction 

In this section the q — 10 Potts model in D = 2, a model which is known to 
show a strong first order phase transition, will be considered; the hamiltonian 
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0.7575 



0.757 



0.7565 



0.756 



0.7555 



0.755 



0.7545 




5e-05 



0.0001 



0.00015 

1/L 2 



0.0002 



0.00025 



0.0003 



Figure 1: (dS(E)/dE) infinite volume extrapolation for the Ising model. 
Fitting function to the data: /(1/L 2 ) = f + f 2 /L 2 . 

of the model is given by 

where the sum is over nearest-neighbours sites and <jj is a Potts spin, taking 
values 1 ... 10. The minimum of the energy is (ferromagnetic ground state) 
and the maximum energy attainable by the system is equal to the number 
of unsatisfied bonds, 2 V. 

The Maxwell construction will be elucidated with the help of figure |5J 
With a standard WLS the DOS of the model for L = 24 in the relevant 
energy range was built. The 1 / -y/ln / criterium [5] has been used to stop the 
updating of \ng(E) at each level of the algorithm, halted when / reached 
the value / min ~ 1 + 10 -7 , and 2 x 10 3 independent runs have been used to 
build the mean. The inverse of the microcanonical equilibrium temperature 
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1.48 




1 37 1 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 1.2 

E/V 

Figure 2: d\og g(E)/dE for D = 2, g = 10 Potts model (L = 24). Green 
band: WLS results. Full blue curve: B-spline fit (12 free parameters) to data 
points. Dotted blue line: the critical value of /? = 1/T as defined by the 
Maxwell construction (see text for the meaning of the labels). 
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at energy E was then estimated by 

0(E) = ^1 = \ng(E + 1) - \ng(E). (12) 

The result is shown as a green band in figure El 

When the two regions A and B have equal areas then the horizontal line 
which defines them at the same time defines the critical value of /3 = 1/T 
at which the first order phase transition occurs. For the two regions to have 
the same area the following condition has to be satisfied: 

JEx aE Je 2 <±E 

This in turns leads to the condition that the double maxima in the energy 
probability distribution, 

P(E) = \ogg(E)-f3E, (14) 

take the same value. E\ and E% identify the locations of the two maxima 
and E 3 is the point of the minimum in between them. It can be shown that 
the common area of the two regions A and B is directly connected to the 
interface free energy between the two phases [121 [T4] . 

In order to perform the integral, a cubic B-spline fit has been used to 
reconstruct the data, obtaining Fwl{L = 24) = 0.10600(31). The error has 
been computed by a standard jackknife procedure, and it has been checked 
that the value obtained is almost completely insensitive to the number of 
free parameters in the B-spline fit; going from 7 to 13 breakpoints (i.e. from 
9 to 15 coefficients) the change in F is very well below the statistical error. 

A subtle point here is the choice of breakpoints locations. The procedure 
used is the one outlined in [13]: instead of choosing equispaced breakpoints, 
the choice of their locations is demanded to a genetic algorithm, and the best 
set of breakpoints is determined by asking for the best overall x 2 . 

The same exercise has then been repeated for several lattice sizes (see 
table [2] for details) up to L = 64. In figure [3] the extrapolation to the infinite 
volume limit is shown. Data have been fitted with a quadratic function: 

F(L)=F + F 1 /L + F 2 /L 2 , (15) 

and the final result, F = 0.0942(8) with x 2 /d.o.f. ~ 0.29, compares quite 
well and is statistically compatible with the analytic result obtained in [14], 
F = 0.094701, and with the multicanonical result of [15J. 
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L 


N 

1 ' run 


F 


24 


2.0 x 10 3 


0.10600(31) 


32 


2.0 x 10 3 


0.10431(22) 


40 


1.8 x 10 3 


0.10301(22) 


48 


1.3 x 10 3 


0.10209(24) 


56 


3.8 x 10 2 


0.10112(35) 


64 


7.5 x 10 2 


0.10022(24) 



Table 2: WLS simulation data for D = 2, q = 10 Potts model. 



0.108 
0.106 
0.104 
0.102 
fe. 0.1 
0.098 
0.096 



0.094 
0.092 

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 

1/L 

Figure 3: Extrapolation to the thermodynamic limit of the interface free 
energy for the D = 2, q = 10 Potts model. Data (red squares) have been 
obtained with a standard WLS and the use of the Maxwell construction. The 
dotted green line represents the analytic result from [14J. 



10 




5 LLR sampling for the Potts model 



The aim of the previous section was only to show that the Maxwell construc- 
tion, together with a good knowledge of the DOS, is a nice tool to numerically 
compute, in the thermodynamic limit, the interface free energy for spin mod- 
els (specifically for the D = 2, q = 10 Potts model), and for this reason the 
WLS computation was not brought to large lattice sizes. In this section it 
will be described the LLRS computation of the same physical quantity The 



L 


AU 


/3 C 


F 


64 


310 


1.421595(8) 


0.11627(19) 


72 


242 


1.422550(9) 


0.11388(22) 


80 


176 


1.423222(9) 


0.11159(24) 


96 


183 


1.424084(7) 


0.10861(25) 


112 


149 


1.424582(6) 


0.10630(25) 


128 


136 


1.424932(6) 


0.10519(26) 


160 


119 


1.425359(5) 


0.10361(29) 


192 


79 


1.425555(9) 


0.10181(50) 


extrap. 




1.426055(9) 


0.09501(54) 



Table 3: LLRS simulation data and results for D = 2, q = 10 Potts model. 

procedure used is the following: the energy range from E min /V ~ 0.28 to 
E max /V ~ 1.08, which is the relevant one for studying the phase transition 
by means of the Maxwell construction, has been considered. The energy 
range has been divided in sub-intervals of size L and for each sub-interval 
a(E) (and whence dlog g(E)/dE) has been computed. The main difference 
respect to the WLS computation regards the fact that Slog g{E)/dE is not 
computed for all values of E but is only sampled, the sampling interval be- 
ing AE = L: the function is fully reconstructed at the end by means of the 
same B-spline fit procedure described above. The overall cost in terms of 
computer time is such that the extension of the computation of F to larger 
lattice sizes is possible in a really reasonable amount of time. The function 
reconstruction can be judged in figure H] (the case L = 128 being taken as an 
example; for all other lattice sizes function reconstruction is of similar qual- 
ity). Once dlog g(E)/dE has been reconstructed, Maxwell construction can 
be used to obtain both the critical value of the temperature at finite volume 
and the value of the interface free energy. Note that the iterating equation 
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1.436 




1 .41 6 1 1 1 1 1 1 1 1 1 1 

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 

ElV 



Figure 4: d log g(e)/dE for L = 128. Continuous line is a cubic B-spline fit 
reconstruction of the data as described in the text. Dotted horizontal line 
is the value of ft c computed via the Maxwell construction. Errors on data 
points are smaller-equal than symbols size. 

for a changes with respect to the Ising case and becomes 

a ^ = a « + 2(6E) 1 +(5Ey {{E * )Kan) - (16) 

Several lattice sizes, from L = 64 up to L = 192, have been considered. 
Simulation parameters (and results for (3 C and F) are in table [3j For each L, 
each run and each value of E, a(E) has been computed by using N t = 100, 
N s = 100 and iV e = 1000. In each case the starting value of a has been taken 
very close to the infinite volume value of /3 C . The first question to address 
is the computation of the critical temperature in the infinite volume limit. 
As can be seen in table [3] and figure [5] is it possible to extrapolate the data 
to the thermodynamic limit as a linear function in 1/L 2 , obtaining a value 
for (3 C = l/T c very close to the analytically known value: j3 c = log(l + vTO). 
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1.427 



1.423 - 
1.422 - 




1.421 - 

1 42 ' 1 1 1 1 1 1 

5e-05 0.0001 0.00015 0.0002 0.00025 

1/L 2 

Figure 5: (3 C extrapolation to the infinite volume limit for the D = 2, q = 10 
Potts model. Data (red squares) have been obtained with LLR sampling and 
the use of the Maxwell construction. Blue points represent data not used 
in the extrapolation. The dotted green line represents the analytic known 
result C = log(l + vTO)- 
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0.12 



0.115 - 



0.11 - 



fe. 0.105 - 



0.1 - 



0.095 - 




0.09 1 1 1 1 1 1 

0.005 0.01 0.015 0.02 

1/L 

Figure 6: Extrapolation to the thermodynamic limit of the interface free 
energy for the D = 2, q = 10 Potts model. Data (red squares) have been 
obtained with LLR sampling and the use of the Maxwell construction. Blue 
points represent data not used in the extrapolation. The dotted green line 
represents the analytic result from [T4] . 

In figure M the infinite volume extrapolation of the interface free energy is 
shown. Note that 1/L finite size effects are much bigger with respect to 
the WLS computation. Nevertheless, discarding the two smaller lattices and 
fitting linearly in 1/L, a good \ 2 around 1 can be obtained, and the final 
result is F = 0.09501(54), which is again compatible with the analytic result 
of [H] and with the multicanonical result of [15J. 

6 Conclusion 

In this work it has been shown that also without a fine tuning of simulation 
parameters (JV e , N t , N s , starting value of a, for example) the algorithm pre- 
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sented in [TJ, which in this paper has been called LLRS, can be successfully 
used with discrete energy models, leading to results for physical quantities 
which in the statistical accuracy reached do not show systematic effects in 
the infinite volume limit. After a short exercise with the Ising model, the 
D = 2 q = 10 Potts model has been considered and the critical value of the 
temperature and of the interface free energy have been computed with a good 
degree of accuracy and a modest computational effort. The main motivation 
for using LLRS instead of WLS, in the case of discrete energy models, is the 
following: 9 log g(E)/dE, which is the key quantity entering in the Maxwell 
construction, is directly computed and can be fully reconstructed for the 
whole relevant energy range by a careful sampling (what is meant here is 
that d log g(E)/dE can be computed for a sample of the discrete energy set 
and then safely interpolated). There is also to mention that in this way one 
can avoid all convergence problems of the energy histogram which have to 
be taken into account when dealing with WLS. 

In this work the possibility to fully reconstruct the DOS, g{E), starting 
from information extracted by LLRS, has not been explored. For continuous 
energy model a simple integration may suffice, and is limited only by dis- 
cretisation effects. For discrete energy model, however, a modest extension 
of the procedure can allow for it. This will be the subject of a future work. 
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